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L  DISCUSSION 

^We  describe  the  SAIC  accomplishments  for  the  final  contract  period. 
Our  previous  annual  report  for  this  contract,  N00014-86-C-2147,  described 
the  work  performed  up  to  March  9,  1988.  Here  the  period  covered  is 
March  10,  1988  to  September  9,  1989.  Because  of  budget  considerations 
the  scale  of  our  program  was  cut  back  and  the  research  performed  was 
undertaken  after  consultations  with  the  COTR  for  the  Laboratory  for 
Computational  Physics  and  Fluid  Dynamics  (LCPFD)  of  the  Naval  Research 
Laboratory  (NRL). 

The  accomplishments  of  the  program  cover  tdiree  major  areas.  They 
include  the  development  of  computer  codes  to  numerically  simulate  (1) 
turbulent  flow  at  a  free-surface  and  (2)  inviscid/viscous  flows  over  complex 
geometries.  The  third  area  concerns  the  development  of  numerical 
algorithms  for  multitarget  tracking.  f  V  j  ( - - 

IL  NUMERICAL  SIMULATION  OF  TURBULENT  FLOW  AT  A  FREE- 

SURFACE 

Introduction 

Much  of  this  year's  effort  was  directed  at  performing  a  direct 
numerical  simulation  of  turbulent  flow  at  a  free-surface.  There  were  several 
distinct  steps  taken  in  this  direction  and  several  distinct  physical  problems 
examined.  All  of  the  work  discussed  below  was  performed  at  NRL  in 
collaboration  with  Drs.  T.  Swean,  R.  Handler,  and  H.  Wang  of  NRL. 

Development  of  Fast  Poisson  Solvers  Using  a  Pseudo-Spectral 

Algorithm 

At  the  center  of  any  Navier-Stokes  solver  using  pseudo-spectral 
methods  are  Poisson  solvers  and  Fast  Fourier  Transforms  (FFT).  For  the 
Fourier  transforms,  the  latest  version  of  the  FFTs  from  Cray  Research  Inc. 
were  used.  For  the  Poisson  solvers,  we  implemented  an  algorithm 
described  by  Orszag  and  Gottlieb  (1977)  and  optimized  it  for  the  Cray  X-MP 
at  NRL.  These  Poisson  solvers  allow  a  single  non-homogeneous  direction 
and  a  single  homogeneous  direction.  Chebechev  polynomials  are  used  as 
expansion  functions  in  the  non-homogeneous  direction  and  trigonometric 
functions  in  the  homogeneous  direction.  In  the  non-homogeneous  direction 
Neumann,  Dirichlet  or  Robin1  boundary  conditions  can  be  used. 


Development  of  a  Two  Dimensional  Navier-Stokes  Solver 

A  two  dimensional  Navier-Stokes  solver  was  developed  as  a  test  bed 
for  the  numerical  algorithm.  The  governing  Navier-Stokes  equations  are 
recast  into  a  4th  order  equation  for  the  vertical  velocity  and  a  2nd  order 
equation  for  the  vertical  vorticity  and  the  continuity  equation  is  solved 
explicitly  in  the  recovery  of  the  streamwise  velocity.  The  equations  are 
numerically  solved  after  they  are  Fourier  transformed  in  the  streamwise  (x) 
and  spanwise  (z)  directions  and  Chebychev  transformed  in  the  vertical 
direction  (y).  For  a  two  dimensional  simulation  the  spanwise  direction  and 
the  equation  for  vertical  vorticity  are  omitted. 

Additionally,  an  Orr-Sommerfeld  solver  was  developed  to  determine 
initial  conditions  for  the  simulations.  Solutions  to  the  Orr-Sommerfeld 
equations  with  the  model  free -surface  boundary  conditions  indicated  that 
the  open  channel,  free-surface  two  dimensional  flow  is  always  stable  in  the 
parameter  range  that  we  considered.  Using  the  Navier-Stokes  solver,  a 
direct  simulation  indicated  that  the  flow  is  also  stable  to  finite  amplitude 
perturbations. 

The  Navier-Stokes  solver  is  also  being  used  to  perform  direct 
numerical  simulations  of  a  vortex  pair  interacting  with  the  free-surface.  In 
this  work,  done  in  collaboration  with  Henry  Wang,  a  velocity  field  vortex  pair 
is  specified  as  an  initial  condition  and  allowed  to  propagate  to  the  free- 
surface.  The  model  free-surface  boundary  condition  was  modified  to  include 
the  effects  of  a  surface  contaminant.  The  addition  of  a  surface  contaminant 
imposed  a  flow  dependent  stress  on  the  free-surface.  This  work  is 
described  in  more  detail  in  Appendix  A,  which  will  appear  in  the  9th 
International  Conference  on  Offshore  Mechanics  and  Arctic  Engineering, 
February  18-25,  1990. 

Development  of  a  Three  Dimensional  Navier  Stokes  Solver 

The  two  dimensional  Navier-Stokes  solver  was  converted  into  a  three 
dimensional  solver  by  the  addition  of  the  equation  for  vertical  vorticity.  A 
strong  emphasis  was  placed  on  both  the  computational  speed,  or  efficiency, 
of  the  computer  code  and  its  versatility  in  terms  of  the  application  of 
boundary  conditions.  The  program  runs  at  40  percent  of  the  theoretical 
speed  of  the  Cray  X-MP.  The  code  was  developed  to  allow  both  slip  and  no 


slip  boundary  conditions  and  either  zero  normal  velocity  or  a  specified 
normal  velocity  (blowing  or  sucking  boundaries  in  the  parlance  of  turbulence 
research).  Much  of  the  effort  of  the  current  year  will  be  directed  towards 
developing  new  boundary  conditions. 

Analysis  of  Fre e-Surface  Turbulent  Flow 

The  incompressible  three-dimensional  Navier-Stokes  equations  are 
solved  for  initial  and  boundary  conditions  approximating  a  turbulent  open- 
channel  flow  of  water  at  Reh  =  3000  where  h  is  the  channel  depth.  Most  of 
the  calculations  were  performed  on  a  32x65x64  grid  in  x,  y,  z  respectively. 
This  resolution  allows  six  grid  points  to  occur  within  the  viscous  sublayer. 
In  the  spanwise  and  streamwise  directions,  the  grid  spacing  is 
approximately  O.lh  and  0.2h  respectively.  These  should  be  compared  to  the 
physically  relevant  scales  of  the  low  speed  streaks  in  the  wall  region:  lh  and 
5h,  respectively  at  this  Reynolds  number.  All  essential  turbulent  scales 
needed  for  the  determination  of  most  statistical  and  flow  structure 
properties  have  been  resolved,  as  can  be  determined  by  comparison  to 
relevant  experimental  and  numerical  results.  A  subgrid  or  large  eddy  model 
has  not  been  used.  The  boundary  conditions  are  periodic  in  all  dependent 
variables  in  the  streamwise  and  span-wise  directions.  No  slip  conditions  are 
used  at  the  channel  bottom  while  the  free  surface  is  approximated  as  a  rigid 
free  slip  surface  with  vanishing  shear. 

A  large  number  of  turbulence  statistics  are  computed  in  the  vicinity  of 
the  free  surface  and  complete  determinations  of  the  balances  of  the  exact 
Reynolds  stress,  turbulence  kinetic  energy,  and  isotropic  dissipation  rate 
equations  are  reported  for  the  first  time.  The  results  show  that  while  the 
turbulence  kinetic  energy  is  preserved  in  the  vicinity  of  the  free  surface,  the 
turbulence  is  redistributed  from  the  vertical  component  into  the  two 
horizontal  components.  The  vertical  vorticity  at  the  free  surface  is 
concentrated  in  regions  elongated  in  the  streamwise  direction.  This 
anisotropic  behavior  leads  to  preferential  redistribution  of  the  turbulence 
kinetic  energy  into  the  span-wise  component  of  kinetic  energy.  The 
balances  of  the  streamwise  and  span-wise  components  of  the  turbulence 
kinetic  energy  reveal  a  reversal  in  sign  of  the  pressure-velocity  correlations 
in  this  region.  A  physical  model  has  been  suggested  to  explain  this  behavior. 
It  is  apparent  from  the  kinetic  energy  and  dissipation  balances  that  there 


3 


exist  two  separate  regions  near  the  free  surface,  a  thin  viscous  layer  and  a 
thicker  zone  wherein  the  redistribution  of  turbulence  is  more  pronounced. 
Near  surface  expansions  of  the  turbulence  kinetic  energy  and  isotropic 
dissipation  rate  are  determined  for  use  in  Reynolds- averaged  turbulence 
models.  A  paper  describing  this  work  is  in  preparation  and  will  be 
submitted  for  publication. 

m.  SIMULATION  OF  INVISCID/VISCOUS  FLOWS  OVER  COMPLEX 

GEOMETRIES 

The  use  of  unstructured  grids  for  the  simulation  of  high-speed  flows 
can  be  found  in  the  literature  (see  references  cited  in  the  appendices).  In 
the  present  research  effort,  we  have  extended  this  technology  to  nearly 
incompressible  flows,  and  applied  the  procedure  to  simulate  inviscid  as  well 
as  viscous  flows  past  submarine  configurations  with  all  its  appendages.  One 
attractive  feature  of  using  triangular  or  tetrahedral  meshes  over  structured 
meshes  is  that  complex  geometries  can  be  easily  represented.  For  example, 
constructing  a  structured  mesh  around  a  submarine  with  all  its  appendages 
will  require  a  tedious  task  of  decomposition  of  the  domain.  In  the  present 
work,  unstructured  grids  are  generated  using  the  advancing  front  algorithm 
of  Lohner.  The  governing  equations  of  flow  are  solved  using  the  finite - 
element  version  of  the  Flux-Corrected  Transport  algorithm  (FEM-FCT). 
Details  of  the  flow  solver  can  be  found  in  the  appendices  referred  to  in  this 
section. 

As  a  first  step.  Euler  and  Navier- Stokes  solutions  were  obtained  for 
axisymmetric  flows.  This  provided  an  excellent  case  to  validate  the 
procedure  employed  and  also  a  base  to  build  models  for  predicting  turbulent 
flows.  The  procedure  was  applied  to  solve  a  model  problem  of  flow  over  a 
sphere  and  the  computed  results  were  found  to  be  in  good  agreement  with 
those  found  in  the  literature  for  both  the  potential  flow  case  and  the  case  of 
viscous  flow  at  Reynolds-number,  Re  =  100.  These  are  part  of  the  paper 
submitted  for  publication  in  the  IJNMFD  journal,  which  is  included  here  as 
Appendix  B.  Having  established  the  correctness  of  the  procedure,  it  was 
then  extended  to  compute  flow  over  the  submarine  hull  configuration.  Grid 
refinement  studies  were  conducted  for  the  inviscid  flow  in  order  to 
establish  the  independence  of  the  flow  solution  to  the  chosen  grid.  Also,  a 
laminar  viscous  flow  solution  for  this  configuration  was  obtained  for  Re  = 


1000.  The  convergence  rate  for  this  problem  deteriorated  considerably,  as 
would  be  expected,  due  to  the  presence  of  the  small  elements  in  the 
boundary  layer  which  are  needed  to  resolve  the  high  gradients  present  in 
the  flow  variables.  Hence,  convergence  acceleration  of  the  numerical 
method  was  investigated  by  appropriately  sub-stepping  the  viscous  diffusion 
terms.  It  was  found  that  this  method  of  convergence  acceleration  did  not 
yield  substantial  gain  because  the  allowable  time-step  for  the  explicit 
scheme  for  low  Mach  numbers  is  limited  by  the  speed  of  sound.  For  explicit 
schemes,  the  allowable  time-step  due  to  the  advective  terms  is  given  by 


Atadv 


< 


A 

Me 


where  A  is  the  minimum  cell  size,  u  is  velocity  and  c  is  the  speed  of  sound. 
For  low  subsonic  flows,  the  allowable  Atadv  therefore  decreases.  Hence,  this 
convergence  acceleration  procedure  should  be  investigated  with  the  barely- 
implicit  correction  (BIC)  scheme. 

The  procedure  was  next  extended  to  solve  three-dimensional  flows. 
Results  were  obtained  for  inviscid  flow  over  the  submarine  with  sail  and 
stem  appendages  at  various  pitch  angles  of  attack.  This  work  was  presented 
at  the  APS  meeting  in  November  1989,  and  an  abstract  of  this  presentation 
is  included  in  this  report  as  Appendix  C.  In  order  to  predict  the  formation 
of  vortices  and  hence  the  ncise  generated  by  them,  it  is  important  to  carry 
out  a  Navier-Stokes  analysis.  Therefore,  the  viscous  diffusion  terms  were 
incorporated  into  the  3-D  version  of  the  flow  solver.  In  the  numerical 
procedure,  these  terms  were  treated  as  a  deferred  correction  in  the  second 
step  of  the  Taylor-Galerkin  procedure.  Preliminary  coarse  grid  results  of 
the  fully  appended  model  at  a  pitch  angle  of  attack  of  10°  show  the  presence 
of  vortices  at  the  junction  of  the  sail  and  the  hull  and  also  at  the  tips  of  the 
stem  planes.  This  configuration  was  also  studied  at  a  yaw  angle  of  attack,  in 
order  to  predict  the  forces  and  moments  that  will  be  involved  in  a 
maneuvering  submarine.  This  effort  will  be  presented  at  the  AIAA  28th 
Aerospace  Sciences  Meeting  and  an  extended  abstract  is  attached  as 
Appendix  D. 
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IV.  ALGORITHM  DEVELOPMENT  FOR  MULTITARGET  TRACKING 


This  section  describes  the  SAIC  effort  in  the  application  of  advanced 
computer  science  algorithms  for  multitarget  tracking.  The  effort  starts  with 
the  notion  that  near  neighbors  finding  algorithms  can  be  used  to  help  track 
objects,  and  it  has  evolved  into  the  identification  and  solution  of  an 
important  problem  that  occurs  in  state-of-the-art  multitarget  tracking 
algorithms. 

The  part  of  the  tracking  problem  we  have  focused  on  is  the 
combinatorial  explosion  that  can  occur  in  gating.  The  previous  year  annual 
report2  gives  an  overview  of  tracking  bottlenecks  and  it  includes  suggestions 
for  efficiency  improvements  in  the  various  parts,  and  in  particular  for  gating. 
A  recent  SAIC  report3  has  detailed  descriptions  and  test  results  of  the 
algorithms  we  have  devised  for  efficient  gating  performance.  The  results  are 
encouraging.  A  copy  of  this  report  is  included  in  Appendix  E.  In  the 
following,  we  outline  the  essential  results. 

We  have  discovered  that  we  can  perform  gating  efficiently  if  we  do 
three  things:  (1)  use  the  characteristics  of  the  Gaussian  correlation  measure 
to  obtain  a  Euclidian  search  radius  from  values  of  the  measure.  This  allows 
the  use  of  efficient  geometrical  computational  methods;  and  it  is  a 
breakthrough  because  the  correlation  measure  is  a  function  of  report  and 
track  covariance  distributions  as  well  as  on  their  position  distributions.  (2) 
Recognize  that  once  you  have  a  search  radius  you  can  use  existing  fast  near- 
neighbors  finding  algorithms  to  pair  the  reports  to  the  tracks.  (3)  In  the 
case  where  the  observation  times  of  the  reports  are  unequal,  use  an  auxiliary 
algorithm  that  makes  "multiple  projections"  of  the  search  data  structures 
and  helps  the  near  neighbors  finding  algorithms  retain  good  overall  scaling. 

Applying  these  techniques  takes  us  one  significant  step  forward  in 
having  a  correlator/tracker  that  can  process  very  large  numbers  of  objects, 
including  SDI  scenario  numbers  of  objects,  in  real  time.  This  is  principally 
so  because  the  algorithms  used  to  do  the  gating  are  designed  for  optimal  or 
near  optimal  scaling,  and  it  can  be  seen  from  analysis  and  data  presented  in3 
that  they  approach  NlnN  scaling. 
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APPENDIX  A 


Direct  Calculation  of  the  Interaction  Between  Subsurface  Vortices  and 
Surface  Contaminant  Distributions  Using  Spectral  Methods 
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Paper  No.  0MAE-90-332 


H.  T.  Wang 

Laboratory  for  Computational  Physics  and  Fluid  Dynamics 
Naval  Research  Laboratory 
Washington,  D.C. 

R.  I.  Leighton 

Science  Applications  International  Corporation 
McLean,  Virginia 


ABSTRACT 

This  paper  presents  a  numerical  calculation  of  the  evolution  of  the 
flow  due  to  a  pair  of  vortices  rising  toward  a  boundary  surface.  A 
spectral  method,  with  transforms  in  Fourier-Chebyshev  space,  is  used  to 
solve  the  two-dimensional  incompressible  Navier-Stokes  equations. 
Three  conditions  are  considered  at  the  surface:  the  fixed  no-slip  and 
shear-free  conditions,  and  a  novel  variable  shear  condition  due  to  the 
presence  of  a  nonuniform  contaminant  distribution.  An  additional  tran¬ 
sport  equation  for  the  contaminant  concentration  is  required  in  this  case. 
The  computation  scheme  is  stable  for  low  and  moderate  contaminant 
concentration  gradient  levels.  Contour  and  line  plots  of  the  main  and 
surface  flows  generally  show  the  expected  result  that  the  variable  shear 
case  lies  between  the  no-slip  and  shear-free  cases. 

INTRODUCTION 

The  use  of  vortex  elements  is  a  powerful  and  convenient  method 
10  model  the  velocity  field  in  a  fluid.  These  elements  may  take  the 
form  of  points,  blobs,  sheets,  or  curved  filaments.  If  bounding  walls 
and  fluid  viscosity  are  neglected,  the  resulting  fluid  velocities  are 
obtained  by  calculating  the  locations  of  the  vortex  elements,  which 
move  at  the  local  velocity,  and  their  induced  velocities.  Thus,  the  prob¬ 
lem  reduces  largely  to  kinematics.  In  the  case  of  a  reflecting  wall,  the 
techniques  of  complex  variables  are  often  used  to  account  for  their  pres¬ 
ence.  In  the  case  of  fluid  viscosity,  it  is  necessary  to  amend  the  calcu¬ 
lation  method  to  account  for  the  diffusion  of  the  vorticity  as  well  as  the 
creation  of  vorticity  due  to  the  no-slip  condition  at  the  wall.  Leonard  [1] 
and  Saffman  and  Baker  [2]  have  written  earlier  surveys  of  the  various 
analytical  and  numerical  vortex  calculation  methods  while  Sarpkaya  [3] 
has  written  a  more  recent  and  comprehensive  survey  of  these  methods. 

Methods  which  consider  the  basic  equations  of  motion  for  the  fluid 
are  less  often  used  to  study  vortex  motions.  These  methods  model 
more  accurately  the  vortex  generation  at  the  bounding  surfaces  and  the 
vortex  diffusion  at  the  expense  of  greater  computational  complexity. 
Examples  of  this  type  of  approach  are  the  studies  of  Ersoy  and  Walker 
[4|  and  Peace  and  Riley  [3],  who  consider  the  case  of  a  pair  of  counter 
rotaung  vortices  approaching  a  bounding  surface.  In  both  studies,  the 
flow  tleld  is  initially  divided  into  an  inviscid  outer  flow  and  a  viscous 
inner  flow.  In  [4],  the  outer  flow  serves  as  a  boundary  condition  for 
the  boundar,  layer  flow  near  a  no-slip  wall,  while  in  [3]  the  outer  flow 


merely  serves  as  the  initial  condition  for  the  solution  of  the  Navier- 
Stokes  equations  for  a  no-slip  as  well  as  a  shear-free  bounding  surface. 
In  both  studies,  finite  differences  are  used  to  obtain  the  spatial  deriva¬ 
tives. 

In  all  of  the  above  studies,  the  principal  interest  is  in  the  vortex 
flow  field  away  from  the  surface,  which  is  taken  to  furnish  fixed  condi¬ 
tions.  In  aerodynamic  applications,  an  important  example  of  such  a 
problem  is  the  flow  field  above  the  runway  due  to  the  vortices  left 
behind  by  departing  or  landing  aircraft.  In  marine  applications,  it  is  of 
interest  to  ascertain  the  effect  of  the  free  surface  on  the  performance  of 
underwater  lifting  surfaces.  More  recently,  improvements  in  remote 
sensing  technology  using  such  techniques  as  synthetic  aperture  radar  and 
infrared  radiometry  make  it  of  interest  to  ascertain  the  wake  features 
around  a  surface  or  submerged  marine  vehicle.  Thus,  the  experimental 
and  theoretical  determination  of  the  surface  elevation  features  of  the 
wake  due  to  submerged  vortices  is  currently  a  field  of  active  interest. 
Examples  of  experimental  investigations  are  those  of  Sarpkaya  [6]  and 
Willmarth,  Tryggvason,  Hirsa,  and  Yu  [7],  who  study  the  surface  wave 
features  caused  by  the  previously  mentioned  case  of  a  pair  of  counter 
rotating  vortices.  Examples  of  numerical  studies  of  this  problem  are 
those  of  Yu  and  Tryggvason  [8]  and  Ohring  and  Lugt  [9],  In  [8|,  a 
potential  flow  combined  vortex/boundary  integral  technique  is  used, 
while  in  [91  the  Navier-Stokes  equations  are  solved  in  curvilinear  coor¬ 
dinates  with  finite  differences  used  to  calculate  the  spatial  derivatives. 
As  may  be  expected,  the  presence  of  the  unknown  free  surface  requires 
an  iterative  procedure  at  each  time  step. 

In  this  paper,  we  present  the  results  of  a  numerical  investigation, 
using  the  Navier-Stokes  equations,  of  the  flow  due  to  a  pair  of  vortices 
approaching  a  free  surface  with  an  unknown  boundary  condition  dif¬ 
ferent  from  the  wave  elevation  case  considered  in  [8,9].  We  use  a  spec¬ 
tral  method,  whereby  the  equations  are  solved  in  Fourier-Chebyshev 
transform  space,  instead  of  the  previously  used  finite-difference 
methods.  It  is  well  known  that  the  Chebyshev  transform  clusters  the 
calculation  points  near  the  wall,  where  the  flow  gradients  tend  to  be 
largest.  In  addition  to  the  fixed  no-slip  and  shear-free  conditions  at  the 
surface,  we  consider  the  novel  condition  of  shear  due  to  the  presence  of 
an  adsorbing  surface  contaminant,  known  also  as  a  surfactant.  In  this 
case,  the  shear  stress  is  variable  and  is  a  function  of  the  constantly 
changing  surfactant  distribution  on  the  surface.  The  modeling  of  this 
condition  requires  a  transport  equation  for  the  surfactant  distribution  at 


the  surface  and  a  coupling  of  the  shear  stress  exerted  by  this  contam¬ 
inant  into  the  boundary  condition  at  the  surface.  Instead  of  using  an 
iterative  procedure  due  to  the  varying  surface  condition,  we  take  advan¬ 
tage  of  our  spectral  approach  to  obtain  a  technique  which  gives  results 
which  are  comparable  in  computer  tune  to  those  for  fixed  conditions. 
This  technique  is.  however,  limited  to  low  and  moderate  contaminant 
shear  stress  cases. 


FREE  SURFACE:  v»0.  u  ■  0,  or  3u/8y  ■  Re  3-rioJ/dx 


We  start  the  body  of  the  paper  by  describing  the  theoretical 
approach.  This  includes  a  description  of  the  Navier-Stokes  equations  in 
rotational  form,  the  derivation  of  the  fourth-order  formulation  which 
implicitly  satisfies  the  troublesome  continuity  condition,  and  the 
transformation  of  the  resulting  formulation  into  Fourier-Chebyshev 
space.  We  give  the  initial  conditions  for  the  vortex  pair,  and  the 
dimensions  of  our  computation  domain.  We  describe  next  the  transport 
equation  for  the  surface  contaminant,  the  computation  technique,  and  its 
limitations.  We  then  present  contour  plots  of  the  velocity  and  vorticity 
distributions  for  the  mam  flow  for  no-slip,  shear-free,  and  contaminant 
surface  conditions.  We  also  present  line  plots  of  the  evolution  of  the 
velocity,  vorticity,  and  contaminant  distributions  on  the  surface  for 
these  same  surface  conditions.  We  conclude  the  paper  by  briefly  sum¬ 
marizing  the  principal  findings. 

THEORETICAL  APPROACH 


0.5  i  0.5 


LOWER  SURFACE:  v-0.  du/dy-0 

Fig.  I  —  Definition  of  Initial  Vortex  Configuration 
and  Computation  Domain 


for  u.  and  using  Eq.(2)  to  write  du/dx  as  -dv/dy,  the  following 
fourth-order  equation  for  v  is  obtained 


±  -  -L 

at  Rt 


By  using  the  initial  vortex  spacing  a.  the  initial  translational  velo¬ 
city  of  the  vortex  pair  Va.  and  the  fluid  density  p.  as  reference  vari¬ 
ables,  the  Navier-Stokes  equations  take  the  following  dimensionless 
form 


where  F  is  the  nonlinear  term  arising  from  u  x  u,  given  by 

F  m  -  _  ^(vu) 

dx2  dydx 


■  w  -L  p2u 

Re 


where  u  is  the  fluid  velocity,  t  is  the  time,  p  is  the  pressure. 
Re  •oaV0ln  is  the  Reynolds  number,  and  p  is  the  fluid  dynamic  viscos¬ 
ity.  For  an  incompressible  fluid,  the  conservation  of  mass  takes  the 
form 

V  u  -  0  (2) 


By  using  vector  identities,  Eq.  (I)  can  be  put  in  the  following  so-called 
rotational  form 


•  vp  +  —  V1 


where  P  •  p  +  u  u/2  is  the  dynamic  pressure  head  and  w-  V  x  u 
is  the  vorticity.  As  noted  by  Hussaini  and  Zang  [10],  the  use  of  this 
form  in  Fourier  collocation  methods,  as  in  our  study,  conserves  kinetic 
energy  and  hence  tends  to  minimize  the  effect  of  nonlinear  instabilities. 

Handler.  Hendricks,  and  Leighton  (11)  point  out  that  a  number  of 
alternate  methods  may  be  used  to  advance  Eqs.  (2)  and  (3)  in  time.  In 
coupled  methods,  the  entire  system  is  considered  at  a  given  time  step. 
In  splitting  methods,  the  time  step  is  split  into  a  momentum  step,  and  a 
step  whereby  the  pressure  is  adjusted  to  satisfy  the  condition  of 
incompressibility,  Eq.  (2).  We  use  here  an  unsplit  scheme,  whereby 
the  troublesome  term  involving  P  is  eliminated  and  the  incompressibility 
condition  is  implicitly  satisfied  by  going  to  a  higher  fourth-order  formu¬ 
lation.  The  approach  implemented  here  is  similar  to  that  proposed  by 
Kim,  Moin,  and  Moser  [12].  For  the  two-dimensional  case  considered 
here,  this  fourth-order  equation  may  be  derived  as  follows.  First,  write 
Eq.  (3)  in  component  form  for  the  velocity  u  in  the  x-direction  and  the 
velocity  v  in  the  y-direction,  where  x  and  y  are  defined  in  Fig.  1 .  Then 
by  taking  d1  /dx1  of  the  equation  for  v  and  -d2 /dydx  of  the  equation 


Numerical  Solution  Procedure 

We  advance  Eq.  1 53  in  time  by  using  the  weighted  implicit 
Crank-Nichoison  method  for  the  linear  term  and  the  weighted  explicit 
Adams- Bash  forth  method  for  the  nonlinear  term,  resulting  in  the  follow¬ 
ing  equation  for  the  value  of  v  at  the  new  n  +  t  time  step  in  terms  of 
values  at  the  previous  n  and  n  -  1  time  steps 


l  — p 

2 Re 


v2v"  +  -J-  (3  r 


where  dU  is  the  size  of  the  time  step.  Since  we  assume  the  flow  not  to 
penetrate  the  upper  and  lower  surfaces  of  our  computation  domain, 
shown  in  Fig.  I ,  the  vertical  velocity  v  is  subject  to  the  following  boun¬ 
dary  conditions 

v  *  0  on  y  »  *  1  (7) 

This  means  that  both  the  upper  and  tower  boundaries  remain  flat.  We 
allow  the  horizontal  velocity  u  and/or  its  derivative  du/dy  to  be  func¬ 
tions  of  x.  Noting  that  du/dx  «  -(dv/dy)  from  continuity,  Eq.  (2). 
the  boundary  conditions  on  u  take  the  following  general  form,  expressed 
in  terms  of  v 

0  »  T~  +  b  #  ~7  “  c  »  on  y  “  *  I  i8> 

dy  dy2 

where  a  ,  and  b  ,  are  given  constants  and  c  .  are.  in  general,  func¬ 
tions  of  x.  Specializing  the  treatment  given  in  (II)  to  the  two- 
dimensional  case,  we  proceed  as  follows  to  solve  for  v'*1.  We 
express  v*  * 1  as  the  following  sum  of  three  partial  solutions 


,/t  - 1 


n  *  1 
P 


(9)  Transformation  to  Founer-Chebvshcv  Space 


t  a  *  v"„* 1  +■  or " 


Each  of  the  individual  solutions  satisfies  the  boundary  condition  given 
by  Eq.  (7),  while  the  two  boundary  conditions  given  by  Eq.  (8)  deter¬ 
mine  the  unknown  coefficients  a*  and  a*.  To  reduce  the  problem  to 
second  order  we  introduce  the  intermediate  variable  f  which  is  related 
to  v  by 


f**1  -  7V*1  (10) 

The  formulation  for  the  particular  solution  v,  which  satisfies  homogene¬ 
ous  boundary  conditions  but  accounts  for  the  nonlinear  term  F  is  given 
by 
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The  formulation  for  the  solution  v  .  which  satisfies  a  nonzero  boundary 
condition  on  y  *  1  is  given  by 
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rr 1  -  o 
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(12b) 
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Finally,  the  formulation  for  the  solution  v.  which  satisfies  a  nonzero 
boundary  condition  on  y  ■  - 1  is  given  by 
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(13b) 
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(13d) 

We  note  here  that  Eqs.  (12)  and  (13)  need  be  solved  only  once  and  give 
O  two  solutions  which  may  be  regarded  as  independent  Green  functions 
whose  coefficients  a*  and  a*  are  adjusted  to  satisfy  the  boundary 
conditions  given  in  Eq.(8). 


In  order  to  solve  the  above  formulation  in  transform  space,  we 
expand  v(x,yj)  given  tn  Eq.(9)  as  a  senes  of  exponential  functions  m 
the  periodic  x  direction  and  Chebyshev  polynomials  in  the  y  direction, 
as  follows 


M/2-1  *  . 

v(x,y.»)  -  £  £  v(m,n,t lexpOlt^r )  7",  (y )  1 14) 

-0  -  -0 

where  M  and  Af  +  1  are  the  number  of  grid  points  in  the  x  and  y  direc¬ 
tions,  respectively,  km  »  2rm  /L,  is  the  mth  wavenumber  in  the  x 
direction,  and  Lt  »  5  (see  Fig.  1)  is  the  length  of  the  computation 
domain  in  the  x  direction.  The  grid  points  xm  and  y„  are  spaced  as  fol¬ 
lows 


mL , 

m  -  0,  1 . ,Vf  —  1 

M 

1 1 3a) 

cos(»n/Af),  n  “0,  1 . Af 

(15b) 

The  Chebyshev  polynomials  are  related  to  the  cosine  functions  by 

T„( y)  ■  cos(/i0),  9  *  cos"'y  (16) 

By  rewriting  the  double  sum  given  in  Eq.  (14)  as  the  following 
single  sum 

M/2-1  „ 

v(x.y.t)  -  £  v(m,y,;)exp(iJr<rr)  (17) 

Ml  "0 


we  reduce  the  problem  jver  the  M  x  (Af  +  1)  physical  points  to  the 
consideration  of  the  .Vf  2  transform  variables  at  each  time  step.  In 
terms  of  the  transform  variables,  Eq.  (Hi),  which  is  the  only  equation 
which  must  be  advanced  in  time,  takes  the  form 


G(y,km.K.K.K")  US) 


In  the  above,  use  has  been  made  of  the  fact  that  the  operation  d/dx  in 
physical  space  corresponds  to  multiplication  by  ii„_in_ transform  space. 
Due  to  the  presence  of  the  nonlinear  terms  Fm  and  Fm  it  is  necessary, 
at  each  time  step,  to  inverse  transform  to  physical  space,  perform  the 
operations  required  to  obtain  F,  and  then  transform  F  to  get  F.  We  per¬ 
form  these  transforms  by  using  standard  Fast  Fourier  Transform  tech¬ 
niques.  Also,  to  avoid  aliasing  errors,  whereby  energy  from  modes 
outside  our  range  of  consideration  is  placed  into  lower  modes,  we  use 
the  well  known  de-aliasing  technique  whereby  we  consider  3M  /2  physi¬ 
cal  points  but  use  only  the  modes  corresponding  to  Af  points. 

We  remark  that  the  transform  variable  f  '  itself  contains  a 
senes  of  N  *  I  Chebyshev  polynomials  (compare  Eqs.  (14)  and  (17)). 
However,  by  using  recursion  relations  which  relate  derivatives  of  a 
Chebyshev  function  of  order  n  to  neighboring  orders,  the  second  deriva¬ 
tive  in  Eq.  (18)  for  the  Af  w  I  grid  points  gives  rise  to  two  ^uasi- 
tndiagonal  matrices  for  the  coefficients  of  the  even  and  odd  Chebyshev 
polynomials  [I3|.  Inversion  of  these  matrices  is  considerably  less  time 
consuming  than  full  matrices  of  the  same  order. 

Boundary  Conditions 

At  the  lower  boundary  y  •  - 1  we  use  the  shear-free  condition 
du  thy  ■  0  to  minimize  tu  effect  on  the  interaction  of  the  vortices  with 
the  upper  boundary,  our  principal  interest.  At  the  upper  boundary 


t 


v  m  l  we  consider  three  conditions.  Two  are  the  standard  fixed  shear- 
tree  and  no-slip  (u  ”  0)  conditions,  os  considered  in  (31.  The  third  is 
the  novel  condition  of  the  presence  of  a  surface  contaminant  which  is 
adsorbed  on  the  water  surface,  or  surfactant.  In  this  case  the  surface 
tension  on  (he  water  surface  y  varies  a  function  of  the  surfactant  con¬ 
centration  <s.  In  the  case  where  <r  varies  with  or  on  the  surface.  Fig.  2 
shows  that  a  shear  stress  occurs,  given  by 

|S.  -  Rt  &&  (19) 


In  terms  of  our  formulation,  which  considers  the  velocity  v  in  transform 
space,  the  above  condition  takes  the  form 


— -  »  Reki,y.  for  m  »  0.  1.2.  .. 

dy 


.  ,  iVf/2  -  l  (20) 


in  Eq.  (22).  Secondly,  we  consider  only  slow  and  moderate  variations 
of  y  with  a.  The  experimental  data  surveyed  by  Skop.  Brown,  and 
Lindsley  [  14)  show  that  the  variation  of  y  with  a  may  be  approximated 
by  several  line  segments,  with  slopes  I A '  \  which  are  typically  less  than 
30  ergs/pg  or  more  than  300  ergs/jig.  We  use  the  absolute  value  sign 
since  A'  is  a  negative  number;  i.e.,  y  decreases  with  increasing  a  We 
find  that  that  we  must  restrict  |  A '  |  to  be  less  than  approximately  60 
ergs/ *ig  to  prevent  significant  numerical  instability.  We  have  been  able 
to  extend  the  range  of  A'  by  using  numerical  damping  techniques  such 
as  incteasing  the  value  of  Rs  or  adding  a  higher  order  diffusion  term 
involving  d*/dx*  in  Eq.(21),  or  (which  is  most  convenient  in  our  spec¬ 
tral  approach)  directly  filtering  out  the  higher  modes  in  Eq.  (22).  We 
have  not  extensively  pursued  these  techniques  since  we  feel  that  exces¬ 
sively  high  damping  will  be  needed  to  stabilize  the  calculations  of  the 
shock-tike  behavor  at  the  high  end  of  |  A '  | .  Also,  as  shown  later,  our 
results  at  the  upper  end  of  the  numerically  stable  |A’|  range  already 
resemble  those  for  the  no-slip  case. 


Eqs.  (19)  and  (20)  indicate  that  it  is  necessary  to  solve  the  following 
transport  equation  for  the  concentration  a  in  order  to  determine  the 
shear  stress  condition  at  the  surface 

da  I  d2a  d(uo)  ... 

dt  “  Rs  4r  ~  dx 

where  0a/d  is  the  Reynolds  number  for  the  surfactant,  and  3  is 
the  surfactant  diffusion  coefficient.  We  solve  Eq.  (21)  implicitly  for  the 
concentration  a  but  couple  it  explicitly  with  the  surface  velocity  u.  By 
using  the  previously  mentioned  Crank-Nicholson  method  for  the  linear 
diffusion  term  and  the  Adams-Bashforth  method  for  the  nonlinear  con¬ 
vection  term.  Eq.  (21)  takes  the  following  form  in  transform  space 


Initial  Conditions 

The  initial  position  of  our  vortices  is  as  shown  in  Fig.  I .  For  the 
Gaussian  vortices  considered  in  our  study,  the  voracity  w,,  i  -  1.2.  of 
each  vortex  is  given  by 


iiiCr.y)  -  — r  exp  -{(a  -  x,)2  +  (y  -  y,)2]/^ 

TT* 


where  X|  »  -0.3,  x*  “  +0.5.  y |  *  yj  ■  -0.75,  and  r  *  0.25  is  a 

measure  of  the  core  size.  The  unit  distance  between  the  vortices  and  a 

value  of  T,  ”  It  give  an  initial  unit  vortex  velocity,  in  accordance  with 
our  nondimensioniization  approach,  discussed  in  connection  with  Eq. 
(1).  The  total  voracity  u  is  given  by 

dv  du  _  2. 

w  =*  - - —  “  2  «i  (24) 

dx  dy 
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It  is  well  known  that  it  is.  in  general,  necessary  to  solve  (he  boundary 
condition  Fq.  '21)  or  (22)  iteratively  with  the  previous  formulation  for 
the  main  flow  until  the  concentration  a  calculated  for  two  successively 
updated  values  of  u  agree  to  within  a  given  error  tolerance  t.  The 
boundary  condition  would,  in  fact,  require  repeated  calculations  of  the 
main  flow,  with  sharply  increased  computer  cost. 


UNIT  DISTANCE 


By  taking  first  derivatives  of  the  above  equation  with  respect  to  x  and.y. 
and  making  use  of  the  continuity  Eq.(2).  we  obtain  two  second-order 
Poisson  equations  for  the  velocities  u  and  v 


Similar  to  what  we  have  described  for  the  time  stepping  procedure,  we 
solve  these  equations  in  Fourier-Chebyshev  transform  space  for  the  fol¬ 
lowing  boundary  conditions.  At  the  lower  boundary  y  *  -1. 
v  •  du  /dy  m  0.  At  (be  upper  boundary  y  m  + 1 ,  v  ■  u  ■  0  for-  the 
no-slip  case,  and  v  ■  du/dy  ■  0  for  the  shear-free  and  surface  con¬ 
taminant  cases.  We  note  here  (hat  the  latter  two  cases  have  identical 
initial  conditions  since  we  take  the  initial  surfactant  distribution  <r(x)  to 
be  uniform. 


4u  . 
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dy 

Fig.  2  —  Shear  Stress  Due  to  Nonuniform  Surfactant  Concentration 


We  have  adopted  a  series  of  approximations  which  render  the  cal¬ 
culation  time  to  be  nearly  the  same  as  for  the  fixed  shear-free  and  no¬ 
il  ip  conditions.  First,  in  order  to  avoid  a  mismatch  between  the  time 
step  requirements  for  the  boundary  condition  and  main  flows,  we  take 
Rs  •  Rt.  We  note  that  numerical  values  for  the  diffusion  coefficient  of 
surfactants  are  not  well  known  and  that  for  the  relatively  high  Rt  used 
in  our  study  (Rt  ■  5000),  the  value  of  Rs  has  only  a  secondary  effect 


NUMERICAL  RESULTS 


Calculation  Parameters 


We  perform  our  calculations  over  a  32  x  33  grid  in  the  x  and  y 
directions,  respectively,  which  corresponds  to  Si  ■  S  -  32.  The 
value  of  Rt  is  taken  to  be  5000.  The  dimensions  of  our  computation 
domain,  in  terms  of  the  initial  vortex  spacing,  is  2  units  in  the  y  direc¬ 
tion.  and  5  units  in  the  x  direction.  These  are  equal  to  those  used  in 
(8].  In  the  following,  we  first  present  contour  plots  of  u.  v,  and  <•>  for 
the  main  flow  and  then  line  plots  of  u,  a,  and  w  at  the  surface  for  vari¬ 
ous  times  t.  We  note  that  r  ■  I  corresponds  to  the  time  required  to 


travel  a  unit  distance  at  die  initial  vortex  velocity  The  results  are 
presented  for  the  no-slip,  shear-free,  and  various  surfactant  cases.  The 
contour  plots  are  shown  for  A'  *  -5.4  ergs/jig,  while  the  line  plots  are 
shown  for  four  values  of  A':  -1.8.  -5.4,  -18.  -54  ergs/ug- 

Contour  Plots  of  Main  Flow 

Figures  3a.  b.  c  respectively  show  plots  of  the  velocity  v  for  the 
shear-free  case  at  t  ■  0.5,  1.5,  2.5.  The  results  are  very  similar  for 
the  no-slip  and  surfactant  cases  with  the  exception  that  the  calculations 
for  the  no-slip  case  had  to  be  stopped  for  /  a  2  due  to  numerical  insta¬ 
bility.  It  seems  that  for  the  high  shear  gradients  in  this  case,  a  finer 
grid  is  necessary.  These  figures  show  the  expected  result  that  at 
t  =  0.5.  the  vortices  are  near  the  lower  boundary,  rise  to  an  intermedi¬ 
ate  position  at  r  »  1.5.  and  flatten  against  the  upper  boundary  at 
r  -  2.5. 
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Fig.  3a  —  Velocity  v  for  Shear-Free  Case  at  i  *  0.5 


The  differences  due  to  the  various  surface  conditions  are  larger  ,n 
the  case  of  u,  and  are  largest  in  the  case  of  u.  Figures  4a.  b.  c  respec¬ 
tively  show  plots  of  u  for  the  shear- free,  surfactant,  and  no-slip  condi¬ 
tions  at  t  »  1.5,  while  Figs.  5a, b  show  plots  of  u  for  the  first  two  cases 
at  t  *  2.0.  Figures  6a,  b,  c  and  Figs.  7a.  b  show  corresponding 
results  for  the  vortictty  u  at  r  »  1.5  and  2.5.  Figures  4  and  6  show  the 
expected  result  that  the  surfactant  case  lies  intermediate  between  the 
shear-free  and  no-slip  conditions.  In  particular,  the  surfactant  case 
tends  to  exhibit  the  main  flow  shape  of  the  shear-free  case  and  (to  some 
extent)  the  surface  shear  characteristics  of  the  no-slip  case.  Figures  5 
and  7  show  that  the  differences  between  the  shear- free  and  surfactant 
cases  are  greater  at  the  later  times.  Now  even  the  main  flow  shapes  of 
the  two  cases  are  perceptibly  different,  particularly  in  the  case  of  j. 

We  also  note  that  the  relative  lack  of  gradients  in  the  shear- free 
flow  case  leads  to  the  computation  being  stable  over  a  longer  time  than 
the  other  two  cases. 
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Fig.  3b  —  Velocity  v  for  Shear-Free  Case  at  t  -  1.5 
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Fig.  3c  —  Velocity  v  for  Shear- Free  Case  at  /  *  2.5 
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Fig.  8c  —  Surface  Surfactant  Concentration  la  at  t  ■  2 
CONCLUSIONS 

We  have  presented  a  spectral  method  involving  transforms  in 
Fourier-Chebyshev  space  to  solve  the  two-dimensional  incompressible 
Navier-Stokes  equations  describing  the  evolution  of  the  flow  due  a  pair 
of  vortices  rising  to  a  surface.  The  use  of  these  transform  methods  con¬ 
verts.  at  each  time  step,  the  calculation  over  the  physical  grid  to  a  series 
of  Poisson  equations  for  the  Fourier  modes.  The  solution  of  each 
Fourier  mode  amounts  to  the  solution  of  sets  of  sparse  algebraic  equa¬ 
tions  for  the  coefficients  of  the  Chebyshev  modes.  By  combining  the 
continuity  equation  with  derivatives  of  the  momentum  equations,  we 
arrive  at  a  fourth-order  formulation  which  eliminates  the  troublesome 
pressure  term. 


At  the  surface  towards  which  the  vortices  are  ruing,  we  consider 
three  boundary  conditions:  the  traditional  fixed  no-slip  and  shear-free 
conditions,  and  a  novel  condition  of  variable  shear  due  to  the  presence 
of  an  adsorbing  surface  contaminant  (or  surfactant),  which  changes  the 
surface  tension  at  the  fluid  surface.  Our  solution  of  the  additional  tran¬ 
sport  equation  modeling  the  surfactant  concentrauon  is  stable,  without 
using  iterative  or  numerical  damping  techniques,  for  low  and  moderate 
magnitudes  of  the  surfactant  surface  tension— concentrauon  slopes. 
Contour  plots  of  the  mam  flow  below  the  surface  show  the  expected 
trend  that  the  surfactant  case  lie*  intermediate  between  the  no-slip  and 
shear-free  cases.  The  difference*  are  quite  small  for  the  verucal  velo¬ 
city  v,  moderate  for  the  horizontal  velocity  u,  and  large  for  the  vorticity 
u.  Line  plot*  of  u,  u,  and  la  (the  deviation  of  the  concentrauon  from 
the  initial  uniform  distribution)  for  the  surface  flow  show  the  manner  in 
which  the  behavior  of  the  surfactant  cases  changes  from  that  of  the 
shear- free  case  to  resemble  the  no-slip  case  with  increasing  magnitudes 
of  the  surface  tension— concentration  slope. 
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Abstract 

This  paper  describes  an  extension  of  previously  devel- 
oped  methodologies  for  solving  the  Euler  and  Navier- 
Stokes  equations  with  unstructured  grids  in  cartesian 
coordinate  systems  [1-5]  to  axisymmetric  coordinate 
systems.  It  is  shown  how  to  arrive  at  a  consistent, 
high-order  formulation  by  proper  choice  of  interpola¬ 
tion  for  the  unknowns.  An  exact  integration  of  all 
integrals  is  performed,  and  the  exact  formulae  are 
derived  and  presented.  Numerical  examples  simu¬ 
lating  both  transient  and  steady-state  flows  in  the 
subsonic,  transonic  and  supersonic  regime  are  given. 
They  demonstrate  the  accuracy  and  wide  range  of  ap¬ 
plicability  of  the  method. 

Introduction 

Axisymmetric  compressible  flow  problems  need  to 
be  simulated  in  many  practical  situations,  including 
flows  in  or  past  bodies  of  revolution  at  sera  angle  of 
attack,  such  as  pipes,  nacelles,  fuselages,  missiles,  as 
well  as  certain  types  of  explosions  and  detonations. 
For  an  axisymmetric  coordinate  system,  the  Navier- 
Stokes  equations  governing  compressible  flows  may  be 
written  as: 

du  ,  dF;  idrn  _st  sf:  i  s. 

ar+ir+;ir  -  T+ir+;ir+T  ■ (la) 

where 


rpv 
rpuv 
rpv 7  +  rp 

rvH 
(16 -d) 
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urr"  +  vrr"  +  rqr  J 

Here  z,r  denote  the  axial  and  radial  coordinates. 
P,  p,  e,  H  denote  the  density,  pressure,  energy  ana  en¬ 
thalpy,  u,  v  denote  the  velocities  in  the  z  and  r  direc¬ 
tion.  Using  Stokes  hypothesis,  the  viscosity  coefficient 
p  and  the  bulk  modulus  A  are  related  by 


and  the  viscous  shear  stresses  and  heat  fluxes  are 
given  by 


„  du  dv 

t".2m5j  >  =  2^  , 

(3a,  b) 

It  n  V  ,,  ( du  dv' 

-<‘U  +  5. 

)  .  (3c,  d) 

,  .or  ,  err 

(3  e,f) 

where  T  and  k  denote  the  temperature  and  thermal 
conductivity  of  the  fluid  respectively.  The  equation 
set  is  completed  by  the  addition  of  the  state  equations 

P  =  (7  -  1)P  [e  -  £  («J  +  v1)] 

r  =  c.[e-i(«’  +  v3)]  , 

t 

(4a,  b) 

which  are  valid  for  a  perfect  gas,  where  7  is  the  ratio 
of  the  specific  heats  and  c,  is  the  specific  heat  at 
constant  volume. 

Multiplication  of  the  system  of  Eqs.  (1)  with  r  yields 


9rU  9rF !  9FI 
9t  +  dz  +  5r 


5  +££l  +  ££L  +  5  (5) 

dz  +  dr 


We  will  denote  the  form  of  the  Euler  equations  as 
given  by  Eq.(l)  as  form  1,  and  the  form  given  by 
Eq.(3)  as  Form  2.  Both  forms  have  been  used  as 


starting  points  for  discrete  approximations.  Form  1 
was  used  by  Kutler,  Chakravarthy  and  Lombard  [6j, 
who  treated  it  as  a  system  of  equations  in  two  dimen* 
sione.  This  straightforward  use  of  form  1  does  not 
produce  a  conservative  difference  scheme,  and  there¬ 
fore  these  authors  employed  a  shock  fitting  scheme  to 
trace  the  shocks.  Form  2  was  employed  by  Deese  and 
Agarwal  [7],  Yu  and  Chen  [8],  and  Woan  [9].  These 
authors  used  this  form  in  Jameson's  two-dimensional 
cell-centered  finite  volume  FL052  code.  Because  the 
scheme  is  cell-centered,  no  problems  appear  at  r  =  0 
(no  nodes  are  placed  there).  However,  problems  are 
expected  at  r  =  0  if  a  node-centered  scheme  is  pre¬ 
ferred. 


Two-Step  Taylor- Galer kin 

The  two-step  Taylor-Galerkin  algorithm  has  been 
used  extensively  for  the  computation  of  both  invis- 
cid  and  viscous  flows  in  two  and  three  dimensions  for 
Cartesian  coordinate  systems  [3-5].  Given  a  system 
of  partial  differential  equations  of  the  form: 


dU  dFj 
9t  +  dxi 


dF* 


(6) 


where  U,  F'  and  5  denote  the  vectors  of  unknowns, 
fluxes  and  source  terms,  we  proceed  as  follows: 

a)  First  sten  :  (Advective  Predictor) 


c- + $f) 


(7) 


b)  Strand  step  ; 


j  W  -£■  2  x  r  dx  dr  ,  rj) 

which  yields  essentially  conservative  form  2.  or 
b)  Take  Conservative  form  2,  interpret  it  as  a  two- 
dimensional  cartesian  problem,  and  incorporate 
it  ‘as  is’  into  an  existing  2-D  code. 

It  is  interesting  to  note  that  whichever  approach  we 
take,  we  always  require  conservative  form  2  in  or¬ 
der  to  obtain  a  consistent,  conservative  scheme.  The 
next  question  that  arises  is  how  to  interpolate  the 
unknowns  involved  in  order  to  obtain  a  discretization 
scheme.  We  can: 

a)  Interpolate  (rp,  rpu,  rpv,  rpe)  by  a  piece- 
wise  linear  approximation.  This  is  the  so-called 
‘group  formulation’.  It  appears  very  economical 
and  simple  to  implement,  but  for  the  limit  as 
r  —  0,  all  derived  quantities,  such  as  the  pres¬ 
sure,  are  not  defined.  They  have  to  be  obtained 
either  using  L’Hopital’s  rule  (which  involves  tak¬ 
ing  derivatives),  or  the  points  lying  on  the  axis 
r  =  0  have  to  be  pushed  to  r  =  c,  where  c  is  a 
small  number.  We  tried  this  option,  but  found 
that  we  always  encountered  numerical  problems 
close  to  the  axis  r  =  0. 

b)  Interpolate  (p,  pu,  pv,  pe)  and  r  by  a 
piecewise  linear  approximation.  This  form  yields 
a  higher  accuracy  in  the  r  -direction  [10]  and  has 
no  problems  at  r  :  0.  The  integrals  that  ap¬ 
pear  in  the  weighted  residual  statement  are  more 
complicated  to  evaluate.  However,  they  may  still 
be  derived  in  closed  form.  For  these  reasons  we 
chose  this  second  form  for  the  spatial  discretiza¬ 
tion  of  the  Euler  equations. 


A U*  =  l/"+l  -  Un  as  At  •  (  S«|"+*  - 

+i-r+f£D  • 


"+i 
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In  both  substepe  the  spatial  discretization  is  per¬ 
formed  via  the  usual  Galerkin  weighted  residual 
method  [3-5].  However,  we  note  that  at  tn+i  = 
tn  +  k  At,  the  quantities  U,  F,  S  are  assumed  as  piece- 
wise  constant  in  the  elements,  whereas  at  t"  ,  tn+l, 
the  quantities  U,  F,S  are  assumed  piecewise  linear. 


Choice  of  Conservative  Form  and  Interpolation 

Having  selected  the  time-marching  algorithm,  we  are 
now  faced  with  the  choice  of  conservative  form.  We 
can  either: 

a)  Take  Conservative  form  I,  and  integrate  consis¬ 
tently,  eg., 


The  First  Step 

Evaluating  all  the  integrals  in  the  weighted  residual 
statement  of  Eq.(5)  and  using  the  notation  defined  in 
Figure  1,  and  the  expressions 


F.i 


r.i  = 


'A+ra  +  rc 
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(10) 


the  following  discretization  for  the  Navier-Stokes 
equations  results: 

Continuity: 
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The  Second  Step 

For  the  second  step,  we  again  evaluate  all  the  inte¬ 
grals  exactly.  Denoting  N'j  as  the  derivative  of  the 
shape  function  N*  with  respect  to  j,  and  Af.  as  the 
consistent  mass  matrix 

Af.  *  J  S*  S*  rdxdr  ,  (15) 

we  obtain  for  the  Euler  equations: 

Continuity: 


Me&gv  =  Af  ^  VOLti  7,1  j  N'a(puvtl  -  rtf) 

.<  L 

+  N*r((jZTi),i-r; f) 
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Energy: 

A/cAge  = 
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Consistent  Mass  Matrices 

A  question  that  arises  from  the  computational  point 
of  view  is  whether  the  consistent  mass  matrix,  which 
is  obtained  by  assembling,  at  element  level,  the  fol¬ 
lowing  exact  element  matrices 


Af.  as 


+  ra 


(20) 


cannot  be  simplified  by  taking  the  average  element  ra¬ 
dius  in  the  integral  ( 13).  This  would  yield  the  element 
matrix 


AfeA£s  Af  ^  VOL. i  7,i  |iV',pti,4  +  S*rpv,i 
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which  is  leas  expensive  to  evaluate.  Our  numeriai  ex¬ 
periments  indicate  that  this  simplification  can  be  em¬ 
ployed  without  loss  of  accuracy.  The  consistent  mass 


matrix  is  solved  iteratively  m  in  the  cartesian  case 
[1*5],  and  again  it  is  found  that  two  to  three  passes 
over  the  elements  are  sufficient  to  raise  the  phase  ac¬ 
curacy  of  the  resulting  scheme  from  second  to  essen¬ 
tially  fourth  order. 

Artificial  Viscosities 

a)  Modified  Lapidus  artificial  viscosity:  The  modified 
Lapidus  artificial  viscosity  [11],  which  proved  success¬ 
ful  for  Cartesian  coordinate  systems,  can  be  extended 
to  the  axisymmetric  case  without  any  further  modi¬ 
fications  by  multiplying  the  element  contributions  by 
their  respective  average  element  radius. 

b)  Maas  diffusion  for  the  FEM-FCT  algorithm:  The 
mass  diffusion  which  is  added  to  the  high-order 
scheme  to  yield  a  monotonic  low-order  scheme  as  part 
of  the  FEM-FCT  algorithm  [5]  can  also  be  extended 
to  the  axisymmetric  case  by  simply  multiplying  the 
element  contributions  by  their  respective  average  ele¬ 
ment  radius. 

Numerical  examples 

A  number  of  numerical  examples  are  given  to  illus¬ 
trate  the  performance  of  the  method  when  simulating 
transient  and  steady-state  problems  in  the  subsonic, 
transonic  and  supersonic  flow  regime.  For  all  the 
steady-etate  problems,  local  timestepping  was  used 
to  accelerate  the  convergence. 

l)  Succragaic  .flow  part  a  iphnc  (rtcufar  rtitcfc  the 
case  under  consideration  corresponds  to  a  free-stream 
Mach  Number  of  Af«o  =  3.0.  For  this  steady-etate 
ease,  only  the  Lapidus  artificial  viscosity  was  em¬ 
ployed  to  stabilise  the  solution.  The  exact  stand-off 
distance  for  the  shock  should  be  of  *  s  1.2 16 A,  where 
R  denotes  the  radius  of  the  sphere  [12].  The  grid  was 
adaptively  remeshed  three  times  [13].  The  final  so¬ 
lution  is  shown  in  Figures  2a-2c.  The  experimental 
stand-off  distance  is  reproduced  exactly  by  the  solu¬ 
tion. 

21  Shock  imoineine  on  a  blunt  body  (transient!:  The 
problem  statement,  as  well  as  the  solutions  obtained 
at  two  different  times  are  shown  in  Figures  3a-3f.  A 
strong  shock  (Atf,  *  10),  coming  from  the  left,  col¬ 
lides  with  the  concave  body  displayed  in  Figure  3a. 
An  adaptive  refinement  scheme  for  transient  problems 
[14]  was  employed  to  resolve  accurately  all  flow  fea¬ 
tures.  The  mash  was  adapted  every  7  times tepe,  and 
two  levels  of  refinement  were  specified.  The  FEM- 
FCT  option  was  invoked  to  maintain  sharp  shock- 
resolution.  The  main  aim  of  this  simulation  was  to 
demonstrate  the  good  phase-accuracy  and  low  nu¬ 
merical  damping  of  the  present  scheme  for  this  class 
of  problems.  As  observed  in  earlier  simulations  of 
this  class  of  problems  [IS- 17]  the  concave  shape  of 


the  body  affects  the  stability  of  the  stand-off  shock 
significantly.  Figures  3g,h  show  the  time-histories  for 
the  pressure  at  a  two  stations  along  the  r  =  0-axis. 
Station  1  (Figure  3g)  lies  at  the  far  right  end  of  the  do¬ 
main,  while  station  7  lies  shortly  behind  the  final  po¬ 
sition  of  the  shock.  One  can  clearly  obeerve  a  damped 
oscillation  for  the  shock  location.  It  takes  many  cycles 
for  the  shock  to  settle  to  its  final  position.  This  be¬ 
haviour,  which  is  not  observed  for  convex  bodies,  was 
also  seen  in  other  numerical  simulations  and  several 
wind-tunnel  experiments  [15-17]. 

3)  Eton  in  aa  U adcru paadcd.  Moalc  (itearix  atattl; 
The  problem  statement,  as  well  as  adapted  mesh 
and  Mach  number  contours  are  shown  in  Figures  4a 
and  4b  respectively.  Several  different  runs  were  per¬ 
formed  for  this  problem.  Some  had  the  FEM-FCT 
option  switched  on,  others  only  employed  the  two- 
step  scheme  described  above.  They  all  showed  the 
existence  of  the  two  shocks  depicted  in  Figure  4b. 
The  run  reproduced  here  was  done  with  a  Lapidus 
artificial  viscosity.  Both  shocks  resulted  from  inad¬ 
equate  nosale  contouring,  as  shown  in  the  contours 
of  the  Mach  number  in  the  region  near  the  throat 
(Fig.  4c).  The  pressure  ratio  across  the  shock  is  sig¬ 
nificantly  lower  than  the  pressure  decrease  through 
the  throat,  though  the  gradients  ace  higher.  During 
convergence  to  steady  state,  the  grid  was  adaptively 
remeahed  three  times.  Tbs  maximum  stretching  ra¬ 
tio  for  the  elements  was  set  to  5  =  9.  A  compari¬ 
son  between  the  measured  and  predicted  radial  dis¬ 
tribution  of  pressure  at  the  exit  plane  is  shown  in 
Fig.  4d.  Significant  scatter  is  shown  in  the  experi¬ 
mental  data,  while  no  data  is  available  in  the  region 
of  the  multiple  shock  system.  Nonetheless,  the  re¬ 
sults  demonstrate  very  good  agreement  over  most  of 
the  exit  plane.  Some  deviation  is  shown  near  the  wall, 
no  doubt  due  to  wall  boundary  layer  effects. 

4)  Flow  past  a  sc  here.  Re=  1QQ  fateadv.  viacQual: 
Steady  viscous  flow  past  a  sphere  at  a  Mach-number 
of  Afa  =  0.1  and  Reynolde-number  of  Re  *  100  pro¬ 
vides  an  important  test  example  to  evaluate  the  ac¬ 
curacy  of  the  present  scheme.  No  artificial  viscosity 
was  added  for  this  subsonic  case.  The  problem  state¬ 
ment,  as  well  as  the  results  obtained,  are  shown  in 
Figure  5.  The  grid  employed  for  this  case,  shown  in 
Fig.  5a,  consists  of  a  structured  portion  divided  into 
triangles  near  the  vicinity  of  the  sphere  and  unstruc¬ 
tured  mash  elsewhere.  From  Fig.  5d,  it  can  be  seen 
that  the  recirculation  tone  extends  upto  1.4D  into  the 
wake,  measured  from  the  center  of  the  sphere.  This 
compares  well  with  experimental  results  [17].  Fig¬ 
ure  5e  shows  the  comparison  of  surface  vorticity  with 
earlier  numerical  results  [18,19],  and  the  agreement 
is  good.  The  flow  separates  at  an  angle  of  approxi¬ 
mately  123  deg. 
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Conclusions 

We  have  described  a  Finite  Element  Solver  for  ax- 
iaymmetric  compressible  flows.  The  Navier-Stokes 
equations  are  advanced  forward  in  time  using  a  two- 
step  Taylor-Galerkin  procedure.  Due  care  was  given 
to  obtain  a  consistent  integration  of  all  variables.  Al¬ 
though  slightly  more  expensive  than  the  equivalent 
2-D  scheme,  the  current  formulation  is  the  only  one 
that  yields  full  second  order  accuracy  for  all  the  un¬ 
knowns  in  both  the  x  and  r-directions.  A  high-order, 
monotonicity  preserving  scheme  is  obtained  by  com¬ 
bining  this  basic  two-step  Taylor-Galerkin  procedure 
with  FEM-FCT  techniques. 

Future  developments  will  center  on  extensions  of 
the  current  explicit  scheme  to  semi-implicit  or  implicit 
schemes. 

Acknowledgements 

This  work  was  partially  funded  by  the  Defense  Nu¬ 
clear  Agency  and  the  Air  Force  Ballistic  Missile  Office 
through  the  Laboratory  for  Computational  Physics 
and  Fluid  Dynamics  of  the  Naval  Research  Labora¬ 
tory. 

References 

[1]  J.  Donea  -  A  Taylor  Galerkin  Method  for  Con¬ 
vective  Transport  Problems;  fnt.  J.  Num.  Math. 
Engng.  20,  101-119  (1984). 

[2]  R.  Lohner,  K.  Morgan  and  O.C.  Zienkiewica  - 
The  Solution  of  Nonlinear  Systems  of  Hyperbolic 
Equations  by  the  Finite  Element  Method;  lot.  J. 
Num.  Metb.  Fluids  4,  1043-1083  (1984). 

[3]  R.  Lohner,  K.  Morgan,  J.  Peraire  and  O.C. 
Zienkiewics  •  Finite  Element  Methods  for  High 
Speed  Flows;  A1AA-85-1S31-CP  (1985). 

[4]  R.  Lohner,  K.  Morgan  and  O.C.  Zienkiewica  • 
An  Adaptive  Finite  Element  Procedure  for  High 
Speed  Flows;  Comp.  Metb.  Appl.  Mecb.  Eng. 
51,  441-465  (1985). 

[5]  R.  Lohner,  K.  Morgan,  J.  Peraire  and  M.  Vah- 
dati  -  Finite  Element  Flux-Corrected  Transport 
(FEM-FCT)  for  the  Euler  and  Navier-Stokes 
Equations;  lot.  J.  Num.  Metb.  Fluid*  7,  1093- 
1109  (1987). 


[6]  P.  Kutier,  S.R.  Chakravarthy  and  C.P  Lombard 
-  Supersonic  Flow  Over  Ablated  Noeetips  Using 
an  Unsteady  Numerical  Procedure;  AIAA  Paper 
78-213  (1978). 

[7]  J.E.  Deese  and  R.K.  Agarwal  -  Calculation  of 
Axiaymmetric  Inlet  Flowfields  Using  the  Euler 
Equations;  AIAA  Paper  83-1853  (1983). 

[8]  N.J.  Yu  and  H.C.  Chen  -  Flow  Simulations 
for  Nacelle-Propeller  Configurations  Using  Euler 
Equations;  AIAA  Paper  84-2143  (1984). 

[9]  C.J.  Woan  -  Euler  Solution  of  Axiaymmetric 
Flows  About  Bodies  of  Revolution  Using  a  Multi¬ 
grid  Method;  AIAA  Paper  85-0017  (1984) 

[10]  P.L.  Roe  -  Error  Estimates  for  Cell-Vertex 
Solvers  of  the  Compressible  Euler  Equations; 
ICASE  Rep.  87-6  (1987). 

[11]  R.  Lohner,  K.  Morgan,  and  J.  Peraire  -  A  Simple 
Extension  to  Multidimensional  Problems  of  the 
Artificial  Viscosity  due  to  Lapidua;  Comm.  Appl. 
Num.  Metb.  1,  141-147  (1985). 

[12]  J.  Zierep  -  Theoretische  Gasdynamik;  G.  Braun, 
Karlsruhe  (1976). 

[13]  R.  Lohner  -  An  Adaptive  Finite  Element  Solver 
for  Transient  Problems  with  Moving  Bodies; 
Comp.  Struct.  30,  303-317  (1988). 

[14]  R.  Lohner  -  An  Adaptive  Finite  Element  Scheme 
for  Transient  Problems  in  CFD;  Comp.  Meth. 
Appl.  Mecb.  Eng.  61,  323-338  (1987). 

[15]  R.  Bastianon  -  Unsteady  Solution  of  the  Flow- 
field  over  Concave  Bodies  -  AIAA  J.  7,  531-533 
(1969). 

[16]  I.O.  Bohachevsky  and  R.N.  Kostoff  -  Supersonic 
Flow  over  Convex  and  Concave  Shapes  with  Ra¬ 
diation  and  Ablation  Effects;  AIAA  J.  10,  1024- 
1031(1972). 

[17]  M.  Van  Dyke  -  An  Album  of  Fluid  Motion  , 
Parabolic  Press;  (1982). 

[18]  R.  Clift,  J.R.  Grace  and  M.E.  Weber  -  Bubbles, 
Drop s  and  Particles  ,  Academic  Press,  (1978). 

[19]  Y.  Rimon  and  S.I.  Cheng  •  Numerical  Solution 
of  a  Uniform  Flow  over  a  Sphere  at  Intermediate 
Reynolds  Numbers;  Phys.  Fluids  12,5,  (1969). 


5 


if  ure  2  :  Supersonic  Flow  Past  a  Sphen 


PRESSURE 


MESH  : 

NELEM=  13116 
NPOIN=  6908 


MACH-NR 


0.00  0.01  0.02  0.03  0.04  0.05  0.00  0.07  0.0ft  0.09  0.10 


J30BIL  * 


RADIAL  DISTANCE.  R/RCW) 


Fig  5.  Radial  Distribution  of  Pressures  at  the  Exit  Plane;  Comparison  between  Predic 
tion  and  Experimental  Data. 
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Figure  3.  Comparison  of  Surface  Vorticity  Distribution  on  a  Sphere. 
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Geometries  Using  a  Finite-Element  Method,* 

R.  Ramamurti,  SAI C  fe  NRL  and  R.  Lohner  GWU  -  The  finite- 
element  method  of  Lohner1  has  been  advanced  to  study  the  flow 
past  complex  3-D  geometries.  In  the  present  investigation,  the 
advancing  front  algorithm3  is  employed  to  generate  the  unstruc¬ 
tured  grids  over  a  complete  submarine  configuration.  A  two-step 
Taylor-Galerkin  procedure  is  used  to  discretize  the  Euler  equations 
of  motion.  The  procedure  was  tested  via  application  to  a  model 
problem  of  inviscid  flow  past  a  sphere  at  Afoo  —  0.2.  Comparison  of 
the  surface  pressure  distribution  with  potential  flow  is  very  good. 
The  procedure  is  then  extended  for  the  simulation  of  3-D  flow  past 
a  submarine  hull  configuration  and  the  results  are  compared  with 
the  axisymmetric  solution.  Flow  past  this  configuration  with  sail 
and  stem  appendages  is  also  investigated  for  various  pitch  angles 
of  attack  to  study  the  asymmetric  flow  properties. 
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Abstract 

A  finite  element  scheme  [1-4]  has  been  advanced  for  solving  the  Euler  and 
Navier- Stokes  equations  with  unstructured  grids  in  both  Cartesian  and  axisym- 
metric  coordinate  systems.  A  two-step  Taylor-Galerkin  procedure  is  employed  to 
discretize  the  governing  equations.  The  accuracy  of  the  scheme  is  validated  by  com¬ 
paring  computed  results  for  flow  over  a  sphere  with  well  known  numerical  results 
and  via  a  grid-  refinement  study  for  an  in  viscid  flow  over  an  axisymmetric  body. 
The  procedure  is  extended  to  solve  three-dimensional  flows  over  submarine  configu¬ 
rations  with  sail  and  stem  appendages.  Convergence  acceleration  for  viscous  flows 
by  sub-stepping  of  the  viscous  terms  is  investigated. 

Introduction 

Numerical  solution  of  flow  past  complex  geometries  is  an  important  tool  for  a 
fluid  dynamidst.  The  use  of  finite  element  methods  using  unstructured  grids  for 
problems  involving  high  speed  flows  can  be  found  in  literature  [1-4].  The  advantage 
of  using  triangular  or  tetrahedral  meshes  over  structured  meshes  is  that  complex 
geometries  can  be  easily  represented.  For  example,  constructing  a  structured  mesh 
around  a  submarine  Mth  all  its  appendages  requires  a  tedious  task  of  decomposition 
of  the  domain.  In  this  work,  unstructured  grids  are  generated  using  the  advancing 
front  algorithm  [5]. 

Most  of  the  conservative  schemes  which  are  extensions  of  1-D  schemes  to  2-D 
and  3-D  through  operator  splitting,  cannot  be  employed  with  unstructured  grids,  as 
the  discretization  stencils  obtained  on  these  grids  are  inherently  multidimensional. 
The  high  resolution  scheme  employed  in  the  present  study  is  based  on  Zalesak’s  [6] 
generalization  of  the  Flux-Corrected  Transport  (FCT)  algorithm  of  Boris  and  Book 


Governing  Equations 

The  equations  governing  the  fluid  flow  are  the  Navier-Stokes  equations,  and  can  be 
written  as  ~ 


4. 4. -  ,•£  .  as  j.  1«SL .  .  A 

dt  dr  dz  ^  r  dx  r*  dr  dz  *  r 

j  =  1,  k  =  0:  axisymmetric  case 
j  =  0,  k  =  1:  3-dimensional  case 

where  U  =  (p,pu,pv,pe)T 

and  Fa  and  Fv  are  inviscid  and  viscous  fluxes  respectively. 

The  equation  of  state  for  an  ideal  polytropic  gas  can  be  written  as 


(1) 


P  =  (7  -  1  )p[e  ~  |  [(/>u)2  +  (/w)2))]  . 


(2) 


For  the  axisymmetric  case  the  system  of  Eqs.  (1)  is  multiplied 


v 


d(rU)  d(rE)  dF  _ 
dt  ^  dx  ^  dr 


d(rEv ) 
dx 


(3) 


Using  this  form  of  the  conservative  equations  can  be  shown  to  be  the  same  as 
integrating  the  system  of  Eqs.  (1)  in  a  consistent  manner.  Moreover,  the  use  of 
conservative  form  represented  by  Eq.  (3)  in  conjuction  with  separate  interpolations 
for  r  and  U  avoids  the  problems  encountered  for  r  =  0  using  a  node-centered  scheme. 


Two-Step  Taylor-Galerkin  Procedure 

The  two-step  Taylor-Galerkin  algorithm  has  been  used  extensively  for  the  com-  • 

putation  of  both  inviscid  and  viscous  flows  in  two  and  three  dimensions  for  Carte¬ 
sian  coordinate  systems  [2-4].  Given  a  system  of  partial  differential  equations  of  the 
form: 


dU  dFi  „  dFt 
dt  +  dx*  ~  Sa  +  dx* 


+  S9  , 


(4) 


where  U,  F*  and  S*  denote  the  vector  of  unknowns,  advective  fluxes  and  advective 
source  terms,  and  F*  and  Sv  denote  viscous  fluxes  and  viscous  source  terms,  we 
proceed  as  follows: 

a)  First  step  (Advective  predictor!: 


-ir+f  .(s.r-fSfr) 


(5) 
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b)  Second  step  : 


A Un  =  U n+1  -Un  =  At-  (Sa|n+*  -  |Q-|n+*  +  ^4|B  +  SJn) 

HT*  f 7T** 


dF' 


dx ' 


5x* 


(6) 


In  both  substeps  the  spatial  discretization  is  performed  via  the  usual  Galerkin 
weighted  residual  method  [2-4].  However,  we  note  that  at  =  tn  +  jAt,  the 
quantities  U,F,S  are  assumed  piecewise  constant,  whereas  at  tn  ,  tn+1,  the  quan¬ 
tities  U,  F,  S  are  assumed  piecewise  linear. 


Convergence  Acceleration 

A  Fourier  stability  analysis  for  the  explicit  scheme  described  above,  shows  that 
the  scheme  is  stable  provided 

C<V1  ~  1  (7) 

where  C  is  the  Courant  number  and  Re&  is  the  minimum  cell  Reynolds  number. 
Convergence  to  steady  state  can  be  accelerated  by  local  timestepping.  Although 
this  local  time-  stepping  strategy  is  efficient  for  inviscid  flows,  convergence  is  rather 
poor  for  viscous  flows.  Hence,  a  sub-stepping  of  the  viscous  terms  is  investigated. 
This  involves  advancing  the  inviscid  fluxes  with  their  maximum  allowable  At  and 
computing  the  corresponding  right  hand  side.  This  inviscid  time-step  is  divided 
into  a  given  number  of  viscous  sub-steps.  The  contribution  from  the  viscous  terms 
to  the  right  hand  side  is  then  computed  and  added  to  the  corresponding  fraction  of 
the  inviscid  right  hand  side.  Complete  details  of  this  procedure  will  be  given  in  the 
final  version  of  the  paper. 


Results 

Axisvmmetric  Flow 

Flow  past  a  Sphere.  Re=  100  (steady,  viscous) 

Steady  viscous  flow  past  a  sphere  at  a  Mach- number  of  M0 0  =  0.1  and 
Reynolds- number  of  Re.  =  100  provides  an  important  test  example  to  evaluate 
the  accuracy  of  the  present  scheme.  No  artificial  viscosity  was  added  for  this  sub¬ 
sonic  case.  The  problem  statement,  as  well  as  the  results  obtained,  are  shown  in 
Figure  1.  The  grid  employed  for  this  case  (Fig.  la)  consists  of  a  structured  por¬ 
tion  divided  into  triangles  in  the  boundary  layer  zone,  and  an  unstructured  mesh 
elsewhere.  From  Fig.  Id,  it  can  be  seen  that  the  recirculation  zone  extends  1.4  di¬ 
ameters  into  the  wake,  measured  from  the  center  of  the  sphere.  This  compares  well 
with  experimental  results  [9].  Figure  le  shows  very  good  agreement  of  computed 
surface  vorticity  with  earlier  numerical  results  [10,11].  The  flow  separates  at  an 
angle  of  approximately  123°. 


Having  established  the  correctness  of  the  procedure,  the  present  scheme  was 
applied  to  solve  flow  past  a  hull-shaped  body  of  revolution.  First,  the  inviscid 
equations  were  solved  on  a  coarse  grid  consisting  of  988  points  and  1807  elements. 
The  results  for  Moo  —  0.2,  are  shown  in  Fig.  2.  In  order  to  establish  the  reliability 
of  the  solutions,  a  grid  refinement  study  is  undertaken.  The  grid  was  refined  using 
the  classic  h-refinement  technique.  The  results  in  terms  of  pressure  contours  and 
velocity  vectors  are  shown  in  Fig.  3.  Figure  4  shows  the  comparison  of  the  surface 
pressure  distribution,  obtained  employing  the  two  grids.  One  can  see  that  the  effect 
of  grid  refinement  is  minimal  on  the  quality  of  the  solution.  This  indicates  that  the 
first  mesh  was  already  quite  adequate. 

Next,  the  procedure  was  applied  to  solve  steady  viscous  flow  past  this  configu¬ 
ration  at  Moo  =  0.1  and  Re  =  1000.  The  grid  employed  for  this  case  consists  of  7276 
nodes  and  14093  elements,  and  is  shown  in  Fig.  5a.  Results  in  terms  of  pressure 
and  vorticity  contours  and  velocity  vectors  are  shown  in  Fig  5b-d.  Correct  trend  in 
surface  pressure  distribution  is  observed.  The  vorticity  contours  show  a  tendency 
for  the  flow  to  separate  in  the  afterbody  region.  This  solution  was  obtained  with 
local  time-stepping  but  without  sub-stepping  of  the  viscous  terms.  Convergence, 
defined  by  reduction  of  residuals  by  three  orders  of  magnitude,  was  achieved  in  5000 
steps.  Convergence  acceleration  for  this  case  is  currently  being  pursued. 


Next,  the  procedure  was  extended  to  3-D  and  an  inviscid  flow  past  a  sphere 
was  chosen  as  the  test  case,  since  axisymmetric  results  from  the  present  study  and 
earlier  results  are  available  for  this  case.  Figure  6a  shows  the  comparison  of  the 
surface  pressure  distribution.  From  this  figure,  it  is  clear  that  the  axisymmetric 
case  compares  very  well  with  the  potential  flow  solution;  the  agreement  of  the  3- 
D  solution  is  fairly  good  except  near  the  two  stagnation  points.  This  discrepancy 
may  be  due  to  the  small  artificial  dissipation  that  was  needed  to  stabilize  the  3-D 
solution  procedure.  Figure  6b  shows  the  pressure  contours  over  the  surface  of  the 
sphere. 


The  procedure  was  extended  to  solve  inviscid  flow  past  a  submarine  with  sail 
and  stem  appendages,  at  a  Mach  number  Moo  =  0.2  and  a  pitch  angle  of  attack 
of  10°.  The  grid  employed  for  this  case  consists  of  410,162  tetrahedra  and  71,524 
nodes  and  is  shown  in  Fig.  7a.  Convergence  to  steady  state  was  achieved  in  800 
iterations,  and  the  results  in  terms  of  surface  pressure  contours  is  shown  in  Fig.  7b. 
Currently,  unsteady  viscous  simulation  of  flow  around  this  configuration  is  being 
pursued,  and  will  be  presented  in  the  final  version  of  the  paper. 
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Figure  i.  Results  for  Flow  past  a  Sphere,  Re  =  100. 

(a)  Grid;  (b)  Pressure  Contours;  (c)  Vorticity  Contours; 
(d)  Velocity  Vectors. 


Figure  i<  Comparison  of  Surface  Vorticity  Distribution  on  a  Sphere. 
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Figure  2  Results  for  Flow  past  a  Hull,  Moo  =  0.2,  Coarse  Mesh, 
(a)  Grid;  b)  Pressure  Contours;  (c)  Velocity  Vectors. 
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Figure  3  Results  for  Flow  past  a  Hull,  M<»  =  0.2,  Fine  Mesh. 

(a)  Grid;  (b)  Pressure  Contours;  (c)  Velocity  Vectors. 


Figure  4  Effect  of  Grid  Refinement  on  the  Surface  Pressure. 


^  Mesh  employed  for  the  viscous  case. 
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I.  Introduction 

Multitarget  tracking  algorithms  have  appeared  in  recent  literature  [1,2,3].  Reference 
1  provided  an  analysis  of  the  scaling  of  times  of  possible  combinatorial  bottlenecks, 
e.g.  gating,  hypothesis  generation,  and  cluster  merging.  Suggestions  were  made  for 
improving  scaling  and/or  the  speed  of  each.  In  this  report  we  focus  on  scaling  and 
speed  improvements  of  gating  algorithms,  whose  combinatorial  problem  we  define  as 
follows:  given  a  set  of  Nr  observations  and  Nt  tracks,  perform  a  matching  algorithm 
such  that  we  identify  all  the  observation- track  pairs,  whose  scores  fall  above  a  certain 
threshold.  For  a  fixed  object  density  or  for  a  variable  density  that  is  sufficiently  low  ,we’ 
verify  that  the  overall  cost  of  making  pairs  will  scale  as  0(Nr\nNt).  For  high  density  and 
observation  reports  with  unequal  timestamps,  we  show  that  the  overall  cost  of  making 
the  pairs  can  be  made  to  scale  better  than  Nt\nNtNr where  (3  is  a  parameter 
known  as  the  search  dimension  and  which  will  be  discussed  below. 

In  the  next  section  we  present  technical  terms  and  definitions.  In  section  III.  we  describe 
the  problem  in  greater  generality  .  We  treat  two  cases:  (1)  the  case  with  either  low 
object  density  or  all  observations  made  at  the  same  time,  and  (2)  the  case  of  high  density 
and  with  observations  made  at  different  times.  In  section  IV.  we  show  that  one  can  use 
near  neighbors  search  algorithms  for  the  first  case,  and  we  discuss  the  derivation  of  a 
search  radius  from  a  pair  score  threshold.  In  sections  V.  and  VI.  we  show  how  to  extend 
these  techniques  to  the  second  case,  with  the  aid  of  an  auxiliary  algorithm.  Scaling 
analysis  is  then  done  for  these  combined  algorithms.  In  sections  VTL  and  VIH.  we  give 
a  brief  discussion  and  conclusions  respectively. 

II.  Preliminary  Details 

By  an  observation  we  mean  a  set  whose  items  are  measured  simultaneously  at  some 
specified  time.  We  call  this  time  the  validity  time  of  the  observation  or  the  ’’time- 
stamp”  of  the  observation.  In  our  discussion  and  in  the  simulations  we  require  that 

’Laboratory  for  Computational  Physics  and  Fluid  Dynamics,  Naval  Research  Laboratory,  Washing¬ 
ton,  D.C.  2037 S 


1 


the  timestamps  all  be  within  a  period  of  time  called  a  "scan”  of  length  Tt.  For  the 
simulation  presented  here,  the  timestamps  are  of  uniform  random  distribution.  A  track 
is  an  estimate  that  in  some  sense  converges  to  the  "true  trajectory”  as  the  number  of 
observations  correctly  matched  to  the  track  increases.  There  are  /?  position  components 
to  each  track  at  any  one  time.  And,  similarly,  among  the  items  of  an  observation  there 
are  (3  position  components.  When  a  geometrical  position  search  is  done,  we  use  all  /? 
components  and  only  these  components.  In  particular,  we  do  not  explicitly  use  the 
observation  validation  time  for  the  searches.  Also,  further  information  may  come  with 
observations,  e.g.  the  location  of  the  measurement  device  and  the  extent  of  its  detection 
volume.  This  additional  information  is  not  used  in  the  techniques  presented  here,  but 
if  used  intelligently  it  could  speed  the  searching. 

We  define  a  score  for  observation-track  pair  (i,j)  as  the  function: 


Sij(dRij,Ti)  <x 


exvjd&Tj'dR) 

detTj1'7 


(1) 


where  Tj  is  the  residual  covariance  of  the  track  j,  and  dR  is  the  vector  position  difference 
of  the  pair,  and  these  arguments  are  taken  to  be  valid  at  the  same  time. 

In  our  simulations  we  use  "near-neighbors-within-a-radius"  algorithms  -  the  acronym  is 
NRA.  The  two  NRA’s  we  use  in  this  study  are:  a  search  tree  [4],  which  in  this  paper 
we  shall  refer  to  as  the  "BLD  tree”;  and  a  modification  of  the  Ordered  Partition  [5] 
which  we  shall  refer  to  as  the  "OP”.  These  are  the  algorithms  of  choice  primarily 
because  of  their  non  statistical  nature,  i.e.  guaranteed  to  find  the  near  neighbors  within 
a  prescribed  radius,  but  also  because  of  their  search  time  scaling. 

III.  Problem  Overview 

In  the  subsequent  discussions  we  assume  that  we  cannot  integrate  the  observations  in 
time.  Given  a  set  of  Nt  tracks  and  a  set  of  Nr  reports,  there  are  at  most  NtNr  scores  5 
which  can  be  formed.  Of  these,  only  a  fraction  q  of  them  will  fall  above  the  threshold, 
and  q  could  be  as  low  as  1  /Nt  or  l/JVr  or  smaller.  Ideally  we  would  only  calculate  the 
qNtNr  scores  and  not  have  any  overhead  for  identifying  these  correct  pairs.  At  worst 
we  would  calculate  the  NtNr  scores  and,  in  addition,  have  significant  overhead  related 
to  this.  An  example  of  a  brute  force  approach  is  the  following  technique:  assume  the 
reports  have  different  times  and  so  for  each  report,  integrate  the  equations  of  motion  of 
the  tracks  to  update  all  tracks  to  the  time  on  the  report  and  to  calculate  Nt  scores.  For 
each  report,  keep  those  scores  that  are  above  the  desired  threshold.  The  dominant  cost 
of  this  is  the  0(NrNt)  score  calculations  and  integrations.  There  are  probably  many 
schemes  for  avoiding  the  brute  force  costs.  The  approach  we  offer  here  is  a  detailed 
outline  of  that  proposed  in  [1]  along  with  data  obtained  from  simulations. 

TV.  Equal  Timestamps  or  Low  Density 

It  is  convenient  to  divide  our  approach  into  two  parts:  first,  how  do  we  use  the  standard 
NRA’s  to  efficiently  do  the  gating  type  of  matching  in  the  straighforward  case  where 
the  sensor  reports  have  equal  timestamps;  and,  second,  how  do  we  use  these  algorithms 


for  the  case  of  unequal  timestamps  without  acquiring  an  0(NtNr )  overhead,  etc.  As  the 
gating  problem  is  posed,  it  would  be  useful  for  the  solution  of  the  first  part  to  calculate  a 
search  radius  from  a  given  score  threshold.  More  specifically,  in  doing  a  search  per  report 
j  on  a  track  database,  we  need  to  determine  the  following:  given  a  threshold,  what  is  the 
radius  squared,  r2(threshold),  such  that  all  the  report-track  pairs  which  are  separated 
by  distances  larger  than  r2( threshold )  will  not  have  scores  above  the  threshold.  Of 
course,  for  a  one  dimensional  gaussian  function  with  a  fixed  T,  the  function  value  is 
uniformly  nonincreasing  as  r2  increases,  and  falls  below  a  threshold  after  r2  gets  larger 
than  some  value,  say  r2(threshold).  It  can  be  shown  that  it  is  possible  to  determine 
analytically  an  r2{threshold)  that  is  independent  of  the  T  distribution  and,  furthermore, 
that  a  useful  one  can  be  found  for  0  >  1  dimensions  [6].  With  such  a  solution  and  its 
obvious  computational  advantages,  all  pairs  that  are  candidate  for  having  scores  within 
a  threshold  are  found,  using  NRA’s,  by  searching  within  a  threshold  radius. 

Though  the  equal  timestamp  behavior  of  the  near  neighbors  algorithm  is  straightforward 
enough,  we  include  Plots  1  and  2  with  data  from  the  simulation.  Plot  1  verifies  the 
0(NtlaNt)  CPU  expense  of  creating  the  data  structure  to  be  searched.  We  have  done 
the  same  simulation  for  two  competing  data  structures,  i.e.  the  OP  and  the  BLD  tree. 
Plot  2  verifies  the  search  time  scaling  for  the  unequal  timestamp  case  with  "relatively 
small”  difference  in  timestamp.  Of  course,  any  difference  in  timestamps  will  force  a 
search  with  a  larger  radius  than  the  case  with  no  timestamp  difference.  Therefore,  we 
consider  this  case  to  be  an  upper  bound  for  the  times  and  scaling  of  searching  and 
scoring  with  equal  timestamps.  It  is  worthwile  to  note  the  nearly  ideal  0(ATrlnJVt) 
search  time  scaling  of  the  BLD  tree,  in  Plot  2,  in  light  of  the  fact  that  Plot  2  is  for 
scenarios  of  fixed  volume  and  both  Nt  and  Nr  were  allowed  to  be  as  as  large  as  32fc. 
That  is,  some  scenarios  were  not  exactly  of  very  low  object  density  abd  the  BLD  tree 
still  performed  very  well. 

V.  High  Density  Unequal  Timestamps:  The  Problem 

When  the  reports  have  unequal  timestamps  what  are  the  best  times  to  which  to  integrate 
the  tracks  for  making  the  potential  pair  matching?  There  are  obvious  answers  if  we  are 
not  concerned  at  all  about  efficiency;  and  also  it  is  obviously  easy  -  from  the  efficiency 
point  of  view  -  if  we  already  know  which  report  to  pair  with  which  track,  but  this 
is  the  problem  we  are  trying  to  solve.  It  is  possible  to  mishandle  this  task  so  that  a 
combinatorial  bottleneck  will  occur  when  matching.  Assume  we  follow  the  textbook 
procedure  of  making  the  tracks  one  dataset  valid  at  the  beginning  or  end  of  the  scan; 
and,  also,  that  we  use  standard  NRAs.  Then  the  r3  argument  to  the  algorithms  is  not 
only  dependent  on  error  bounds  from  a  score  threshold,  an  r2(threahold)  component, 
but  also  on  bounds  on  positions  determined  by  the  location  possibilities  of  the  objects 
due  to  their  dynamics  and  to  the  time  differences.  The  latter  factor  can  be  orders 
of  magnitude  larger  than  the  first.  It  is  not  difficult  to  imagine  situations  in  which 
dynamics  and  object  density  will  cause  the  near  neighbors  finding  algorithms  to  return 
much  of  the  data  in  the  track  dataset  so  that  again  we  approach  the  0(NtNr)  scaling. 
And  this  horrendous  scaling  can  occur  not  only  for  the  score  calculations,  but  also,  in 
the  neighbors  finding  section.  This  is  so  because  the  time  of  the  NRAs’  execution  is 
dependent  on  the  number  of  neighbors  it  returns,  which  in  turn  is  dependent  on  r3,  as 
will  be  argued  later.  In  Plot  4  we  have  redone  Plot  2  except  that  the  scan  length  is 


increased  by  a  factor  of  five  and  therefore  the  average  search  radius  is  also  increased. 
In  this  case  neither  of  the  two  search  algorithms  scale  as  0(Nr\nNt).  The  search  times 
increased  by  a  factor  of  about  2  and  5,  for  the  OP  and  the  BLD  tree  respectively,  and 
the  number  of  near  neighbors  returned  was  increased  by  a  factor  of  about  10. 

VI.  High  Density  and  Unequal  Timestamps  :  Solution  Description  and  Analysis. 


Part  of  the  motivation  for  our  approach  lies  in  the  following  observations:  (1)  that 
some  data  structures  for  NRAs  are  cheap  to  make  -  CPU  wise,  and  (2)  that  by  making 
various  copies  of  the  track  data  structure  -  the  various  copies  valid  at  different  times  - 
a  tradeoff  can  be  made  between  the  time  spent  on  the  creation  of  the  data  structures 
and  the  time  spent  on  searching  and  scoring.  The  payoff  from  spending  more  time  in 
creating  data  structures  increases  as  the  number  and  density  increase  and  therefore  may 
reduce  the  scaling  of  overall  gating.  Consider  the  case  of  a  luge  number  of  objects  with 
extremely  high  density  and  a  scan  length  Tt  so  that  the  r3  argument  to  the  NRA  gives  a 
search  volume  comparable  in  size  to  the  extent  of  the  object  distribution.  Also,  suppose 
the  object  distribution  is  equally  dense  and  random  in  each  direction.  Then,  naturally, 
many  if  not  all  objects  will  be  returned  as  candidates  by  the  NBA  and  we  are  near  the 
0(NtNr)  scaling  that  we  want  to  avoid. 


If  we  make  Ni  track  data  structures  with  NRAs  valid  at  A,  equally  spaced  time  intervals 
within  the  scan  of  length  T„  then  any  report  would  be  at  most  T,/{2xNi)  time  units 
away  from  a  track  data  structure  (TDS).  The  average  radius  to  the  NRA  is  decreased 
by  the  factor  Ni  -  compared  to  the  case  where  we  have  one  copy  at  the  middle  of  the 
scantime  -  and  the  volume  extent  as  well  as  the  average  number  of  candidates  returned 
is  smaller  by  (1  /Ni)3  in  the  isotropic  dense  limit  3D  case.  That  is,  in  this  limit  the  total 
number  of  N  objects  in  Nr  spheres  each  of  radius  r,-  =  %VmaxT,/Ny3  is  given  by 


A  Nr/2Nj  jr  rp 

N  =  A*  £  (2) 

J  .x.1  Nt 

where  p  is  the  track  density  and  Vmax  is  the  speed.  The  sum  has  an  exact  solution  as 
a  function  of  Nr,  but  in  the  large  Nr  limit  we  have 

N  =  />(4x/3)21Vi(Vm„rf/lVt1/3)3(3~-)(^-)3+1  (3) 

and  for  a  fixed  Nt,  Nr  and  T,  then: 


N  «  (1  /Ni)3  (4) 

Overall  CPU  time  could  be  reduced  since  scoring  time  is  directly  proportional  to  the 
number  of  near  neighbors  found  as  is  searching  time  -  though  exactly  how  depends  on 
the  particular  NRA.  Though  our  simulation  is  somewhat  degenerate  in  one  of  the  3 
space  dimensions,  for  only  32 K  objects  and  T,  =  10s,  the  number  of  near  neighbors 
returned  decreased  by  an  average  factor  of  9.8  when  Ni  was  changed  from  1  to  5.  Also, 
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for  Ni  =  5,  the  number  of  near  neighbors  returned  was  decreased  by  a  factor  of  62.7 
compared  to  the  case  of  1  TDS  placed  at  the  beginning  of  the  scan.  Plot  5  shows  how 
the  scaling  of  the  search  time  behaves  for  the  same  scenario  as  in  Plot  4,  except  that 
now  Ni  —  5  for  the  datasets  of  4 K  and  32 K  tracks.  The  best  fit  lines  approach  the 
N\nN  scaling  line  more  closely.  Also,  for  the  same  runs,  Plot  5  shows  that  with  multiple 
copies  the  scaling  for  scoring  was  less  than  0(N\rx.N),  and  the  data  shows  it  was  almost 
linear.  Of  course,  a  price  was  paid  in  creation  of  the  multiple  copies  and  data  structures. 

Was  the  price  paid  for  the  multiple  data  copies  too  high  so  that  the  overall  cost  would 
still  scale  worse  than  IVlnJV  ?  Let  k  be  the  average  number  of  near  neighbors  returned. 
Then  the  total  CPU  cost  can  be  modeled  as 


T(NU  Nr,  Ni , k, )  =  CtNiNt  +  C2NtNt\nNt  +  (C^NAnNt  +  C3bNrk)  +  C4kNr  (5) 

where  the  terms  on  the  right  hand  side  of  equation  (5)  give  respectively,  the  cost  for 
integrating  the  tracks  to  the  desired  time  of  the  data  structures,  the  cost  of  making  the 
data  structures  (BLD  tree  case),  the  cost  of  searching  the  appropriate  data  structure  for 
each  report(BLD  tree  case),  and  the  cost  of  scoring  the  pairs.  Of  course,  k  is  a  function 
of  the  density  and  the  total  volume  searched: 

k  =  k(p,r2).  (6) 

where  p  =  p(Nt,Vol )  and  r2  =  r2(T„  Ni,T ,  Vmax).  In  the  case  of  constant  volume  -  the 
case  we  actually  modeled  -  a  fixed  T$,  dense  and  isotropic  and  in  which  T  provides  a 
correction  to  r  of  lower  order  than  that  determined  by  the  dynamics,  then 


k<xk(NuNi)<xNt/{N?).  (7) 

Now  we  can  rewrite  equation  (5)  as 

T(Nr,Nt,Ni)  =  CiNNt  +  CtNiNt]nNt+ 

C^NMNt  +  C&Nr&)  +  CA,NrA).  (8) 

The  function  above  has  a  minimum,  which  can  be  verified  to  be 


Nimtn  =  Wr(Cw  +  Cv)/(Ci  +  CalnW,))*1^1”  (9) 

and  substituting  equation  (9)  into  (8)  we  find  an  upper  bound  on  the  leading  scaling 
term  goes  as  0(NtlnNtN^l^l3+1^).  Of  course,  with  random  timestamps  within  a  scan, 
only  in  the  assumed  dense  limit  will  it  scale  this  poorly  as  compared  to  the  NrlnNt 
scaling  of  the  low  density  case. 


VII.  Discussion 


Not  only  is  the  value  of  JV,-,  for  a  given  Nt  and  Nr,  software  dependent,  but  it  is  also 
dependent  on  the  computer  used.  For  example  a  vector  computer  with  good  floating 
point  hardware  will  no  doubt  reduce  C v  over  a  generic  computer,  but  because  C\  can  be 
reduced  by  an  even  larger  factor  since  vectorization  of  the  integration  of  the  equations  of 
motion  will  be  straightforward  and  possibly  complete.  Thus  the  values  of  the  constants 
will  need  to  be  "in  house  empirical”  in  every  hardware/software  setup. 

Because  only  scaling  data  was  provided,  we  should  mention  that  the  code  currently  takes 
about  1  minute  on  a  Sun  3/260  to  process  a  case  in  which  Nt  =  Nr  =  4 K.  By  processing 
we  mean  the  time  to  generate  all  the  data  and  to  find  and  score  the  candidates  twice  • 
once  for  each  NRA. 

Vni.  Conclusion 

In  this  paper  we  have  presented  an  efficient  approach  to  gating  in  multitarget  tracking.  A 
technique  was  presented  to  (1)  handle  the  case  of  fixed  object  density,  or  sufficiently  low 
density,  with  equal  or  nearly  equal  observation  timestamps,  and  (2)  to  handle  the  case 
of  high  density  with  observation  reports  of  unequal  timestamps.  For  the  former  case  we 
verified  that  the  overall  cost  of  making  observation-track  pairs  scale  as  0(NrlnNt).  For 
the  latter  case  we  showed  that  the  overall  cost  of  making  the  pairs  can  be  made  to  scale 
better  than  NtlnNtN^1^l3+l^\  In  both  cases  we  used  “near-neighbors-within-a-radius” 
search  algorithms,  and  in  the  latter  case  we  introduced  a  new  auxiliary  algorithm  to 
achieve  the  stated  gating. 
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Search  Time  Scaling  of  Data  Structures.  Search  for  neighbors 
of  N  reports  among  data  structures  with  N  tracks;  return  Pairs. 
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Search  time  scaling. 
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